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Abstract. In open double-well Bose-Einstein condensate systems which balance 
in- and outfluxes of atoms and which are effectively described by a non-hermitian 
■PT-symmetric Hamiltonian PT-symmetric states have been shown to exist. PT- 
symmetric states obey parity and time reversal symmetry. We tackle the question of 
how the in- and outfluxes can be realized and introduce a hermitian system in which 
two PT-symmetric subsystems are embedded. This system no longer requires an in- 
and outcoupling to and from the environment. We show that the subsystems still have 
PT-symmetric states. In addition we examine what degree of detail is necessary to 
correctly model the PT-symmetric properties and the bifurcation structure of such 
a system. We examine a four-mode matrix model and a system described by the 
full Gross-Pitaevskii equation in one dimension. We see that a simple matrix model 
correctly describes the qualitative properties of the system. For sufficiently isolated 
wells there is also quantitative agreement with the more advanced system descriptions. 
We also investigate which properties the wave functions of a system must fulfil to 
allow for PT-symmetric states. In particular the requirements for the phase difference 
between different parts of the system are examined. 


PACS numbers: 03.75.Kk, ll.30.Er, 03.65.Ge 

1. Introduction 

In conventional quantum mechanics hermitian operators are used to describe closed 
quantum systems. These operators allow only for real eigenvalues, which can represent 
physical observables. Since systems in the real world are hardly ever completely isolated, 
the environment must be taken into account. Due to a lack of knowledge about the actual 
layout of the environment of a system or because the environment is too complicated 
to be completely calculated, one can effectively describe such systems as open quantum 
systems as long as the interaction to the environment is known. Such Hamiltonians 
are often no longer hermitian. The interaction with the environment, e.g. gain and loss 
of the probability amplitude, that is the wave function, can be expressed by complex 
potentials [^. These Hamiltonians in general do not have a real eigenvalue spectrum. 
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A special class of non-hermitian operators was investigated by Bender and 
Boettcher in 1998 [^. For certain parameter ranges these operators also had purely 
real eigenvalue spectra. The origin of the special property can be traced back to the 
■PT-symmetry of the operator, where the PT-operator consists of the parity operator P 
and the time reversal operator P. The parity operator exchanges x —)■ —x and p —)■ —p. 
The time reversal operator replaces p —)■ —p and i —)■ —i. A PT-symmetric system has 
a Hamiltonian which fulhls [H,VT] = 0. For a system with 

H = -A + V (1) 

the position space representation of the potential must obey the condition 

V{x) = V*{-x), (2) 


i.e. the real part of the potential must be an even function in the spatial coordinate 
and the imaginary part must be an odd function. PT-symmetric systems have been 
studied theoretically for quantum systems j3]-[^. However, the concept of PT-symmetry 
is not restricted to quantum mechanics. Indeed, the experimental breakthrough was 
achieved in optical wave guides by Riiter et al when in such a system the effects 
of PT-symmetry and PT-symmetry breaking were observed. This has led to a still 
increasing interest in the topic (8|-[TT|, and PT-symmetric systems have also been studied 


in microwave cavities 12 , electronic devices 13 14 , and in optical 15 -24 systems. Also 
in quantum mechanics the stationary Schrodinger equation was solved for scattering 
solutions El and bound states IS]. Note that it was shown in [25l that the characteristic 


PT-symmetric properties are still found when a many-particle description is used. 

it was suggested that PT-symmetry could also be realized in quantum 


In 26 


systems, namely in Bose-Einstein condensates (BECs). The BEC would be captured in 
a symmetric double-well potential where particles are gained in one well and lost in the 
other. This loss and gain can then be described by a complex potential coupling the 
system to the environment. 

The time-independent solutions (see Appendix A) of such a PT-symmetric double¬ 
well system can in the simplest possible case (2^ be described by the matrix 


-glA? -h 


-5'|^2p + i7 



= h 



( 3 ) 


where ‘ipi and ■02 represent the occupations of the two wells with atoms in the condensed 
phase and fi is the chemical potential. This description can be derived from a non- 
hermitian representation of a many-particle Bose-Hubbard dimer 28 . The off-diagonal 


elements v of the matrix describe the couplings between the wave functions in the 
two potential wells. The diagonal contains a nonlinear entry introducing the particle- 
particle interaction described by an s-wave scattering process. Its strength can be 
changed via the parameter g which is proportional to the s-wave scattering length and 
its physical variation can be achieved close to Feshbach resonances. In comparison to 
the original model from 27 the replacement g —)■ —g is introduced to be consistent with 
the other models which will be shown later on. In addition the diagonal contains an 
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Figure 1. Analytic solutions for the chemical potential Q of the two-dimensional 
matrix model described in (|^. The coupling strength v = 1, and the nonlinearities 
g = 0 in a), g = 1.4 in b) and g = 2.6 in c) are used. The analytically continued 
solutions are plotted using dashed lines. 


imaginary term with the parameter 7 . This term models a particle gain in one well and 
a particle loss in the other. This gain and loss is provided by the (not further described) 
environment. The wave functions consist of two complex values and contain no spatial 
information. Therefore the parity operator V, which normally exchanges x with —x, 
exchanges "01 with 'ip 2 and vice versa. It is also assumed that the potential wells are 
isolated enough such that the nonlinear interaction between and 'ip 2 can be neglected. 

The system (|^ is solved analytically 27 for wave function vectors ip which are 
normalized to one. The chemical potential reads 


Rs 

Ra 


■| ± \/u2 -72, 


-^±71 


4y2 


g2 _|_ 4^2 


- 1 . 


(4) 


The values in (|^ are the "PT-symmetric solutions, and the PT-broken solutions 
of the system are denoted All solutions are shown in hgure For small 7 the 
system without nonlinearity {g = 0) shows only PP-symmetric states with real chemical 
potential /i G M as can be observed in hgure [^. These states pass through a tangent 
bifurcation at 7 = 7 c = 1, and two PT-broken states emerge. For 7 > 7c only PT- 
broken states with a complex chemical potential fi E C exist. 

For a nonlinearity g > 0 the bifurcation in which the two PT-broken states are 
created moves to a smaller value of 7 on one of the PT-symmetric branches (compare 
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figure &)■ A pitchfork bifurcation is formed. Thus for nonzero values of g there is 
an additional parameter region for 7 , in which "PT-symmetric and PT-broken states 
exist simultaneously. When the nonlinearity is increased even further {g > 2) we see in 
figure that the pitchfork bifurcation is no longer present and the PT-broken states 
exist for all values of 7 . A thorough examination of the bifurcation structure and of the 


associated exceptional points can be found in 29 


The matrix model does not take the spatial extension of the system into account. In 
general BECs can be described by the nonlinear Gross-Pitaevskii equation 30 . Often 


5 functions have been used to gain a deeper insight |4|[^ [3l| - |4l] . Therefore a simple 
model to include spatial effects describes the potential with double-^ functions 42 . In 


this system two h-wells exist at the positions x = ± 6 . While both of these wells have 
the same real depth they possess antisymmetric imaginary parts. That is, one well has 
a particle gain and the other has an equally strong particle drain. The potential fulfils 
the PT-symmetry condition (|^ and the corresponding Gross-Pitaevskii equation is 

— ^|J"{x) — [(1 - 1 - i'y)6{x -\-b) -+■ {1 — i'y)5{x — b)] ipix) — g\'ip{x)\^il){x) = iJiil){x). (5) 

In this system PP-symmetric solutions and PT-symmetry breaking were found. 


In 43,44 a similar two well system was examined in much greater detail by using 


a more realistic potential well shape. The Gross-Pitaevskii equation of such a BEG can 
be written as 


(-A + V{x)- g\i){xR)\^)i) = 


with the complex potential 


V{x) = -x^ + To 




with p = 


a 


( 6 ) 


(7) 


21n(4To^cr) 

containing the BEG in a harmonic trap divided by a Gaussian potential barrier into two 
wells. The parameter p is chosen in such a way that the maximal coupling between the 
subsystems occurs at the minima of the potential wells. The stationary states show the 
same general behaviour as those in the matrix model. 

All descriptions so far used complex potentials to effectively describe the 
environment. Therefore only the PT-symmetric part of the whole system was described 
in detail while the concrete layout of the environment itself was not specified. We will 
now discuss how it might be possible to embed such a PT-symmetric two-well system 
into a larger hermitian system and therefore explicitly include the environment into our 
description. 

where 


As a first step in this direction a hermitian four well model was used 45,46 


the double-well with in- and outgoing particle fluxes is achieved by embedding it into 
the larger system. The two outer wells have time-dependent adjustable parameters 
namely the potential depth and the coupling strength to the inner wells. By lowering 
and raising these wells a particle gain and loss in the two inner wells can be obtained, 
which exactly corresponds to the loss and gain in the non-hermitian two-well model. 
However, the PT-symmetric subsystem of the inner wells loses its properties when the 
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well which provides the particle gain is depleted. A second possible realization was 
suggested in 47 , where the wave function of a double-well potential was coupled to 


additional unbound wave functions (e.g. one ingoing and one outgoing) connecting the 
gain and loss of the system with a reservoir. These auxiliary wave functions replace the 
previously unknown environment of the system. 

In this paper we propose an additional way of realizing a PT-symmetric two-well 
system by extending the approach used in 47 . We couple two stationary bound wave 


functions, where each of them has the shape of that of the corresponding "PT-symmetric 
system and their combination results in a hermitian system. The influx from one system 
originates from the second and vice versa. By tuning the coupling strength between the 
two systems we will be able to control the gain and loss in the subsystems. In contrast 
to 47 our systems are closed and do not require incoming or outgoing wave functions 
or time-dependent potentials. We will show that for suitable states the subsystems are 
indeed PT-symmetric, however, also PP-symmetry breaking can be observed. 

In section |2] a four-dimensional matrix model will be constructed similar to the 
model ([^ to observe the general structure of the eigenstates and to determine their 
PP-symmetric properties. For this model analytical solutions can be found. In a next 
step a Hamiltonian is constructed to combine two subsystems with a spatial resolution 
in one dimension for the wave function similar to the double-Ppotential used in ([^. In 
these systems effects which depend on the shape of the wave functions can be observed. 
We will examine which detail of description is necessary to capture the PP-symmetric 
properties of the system and the bifurcation structure. Since a model with double- 
Ppotentials is only a rough approximation of the reality we will also introduce an 
additional system. This system is constructed by coupling two subsystems of the form 
(§. It not only has an expanded wave function, which resolves spatial information, but 
also possesses more realistic extended potential wells. In addition the coupling between 
the two modes takes place over an extended area of space and is not confined to the 
locations of the Pwells. 

Subsequently we will compare the results obtained with the different descriptions in 
section]^ We will also compare the bifurcation structure with the model ([^. In addition 
the influence of the phase difference between the two subsystems on the stationary states 
will be determined. A summary and discussion of the results is given in section 


2. Coupling of two two-well potentials in one hermitian system 

In hgurej^the layout of two coupled two-well systems is sketched. The two subsystems 
are labelled A and B and each contains two wells with the labels 1 and 2. In the drawing 
the potentials of the wells are extended. This corresponds to an ansatz as shown in (|^ 
and ([^ and will be one of the systems studied in this work. Each of the wells is coupled 
to its counterpart in the other subsystem. The coupling strength is described by the 
parameter 7 . Since the strength of the in- and outcoupling is also determined by the 
wave function of the other subsystem, "PT-symmetry can only exist for both subsystems. 
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Subsystem A 


Subsystem B 


1 2 




M 


-i? 



1 2 


Figure 2. This sketch illustrates how two double-well subsystems are combined into 
a closed hermitian system. The coupling and description of the wells is given with a 
varying degree of detail for the different systems discussed in this paper. 


There is no "PT-symmetry for arbitrary states but only for states with an appropriate 
symmetry. As mentioned in the introduction we will investigate the setup in three 
different degrees of detail. 

There exist various other systems which have four distinguished modes. A family 
of such systems named plaquettes was examined 48,49 . These systems are seen as a 


first step towards building PT-symmetric lattice systems 50,51 . The plaquettes exist 


in various conhgurations which differ in the coupling between the sites. In contrast to 
the model proposed in this paper the gain and loss in these plaquettes is still provided 
by non-hermitian terms. 


2.1. Four-dimensional matrix model 

In a first step we construct the four-dimensional hermitian matrix model. Therefore we 
place two matrices of the shape ([^ on the main diagonal blocks in our new matrix M 
and remove the terms which couple the system to the environment. They are replaced 
with coupling terms in the off-diagonal 2x2 matrix-blocks, i.e. 


M = 



V 

-iy 

0 

\ 

V 

-5'|^A,2p 

0 

+iy 


+i7 

0 


V 


0 

-iy 

V 

-5'IV’B,2p 

/ 


( 8 ) 


with the wave function 




(9) 


The elements of 'ip are four complex values with no information about the spatial 
extension of the wave function. The first two values 'ipA,i,i’A ,2 G C represent the 
wave function amplitudes in the double-well potential of subsystem A while the values 
V’n,i)'0B,2 £ C represent the amplitudes in the subsystem B. Therefore in this context 
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the parity operator V exchanges -^a,! with i/ja ,2 and with V'Bg- The two diagonal 
submatrices will form our subsystems A and B each with two wells indicated by the 
indices 1 and 2. The hrst well of subsystem A is coupled via = —iy to the hrst well 

in subsystem B. The hrst well of subsystem B is coupled via = iy to therefore 

keeping the matrix hermitian. The second wells are coupled in a similar manner but 
with opposite signs. Note that the coupling terms do not yet guarantee a symmetric 
gain and loss in a subsystem since the gain and loss depend also on the value of the 
wave function of the other mode. 

The coupling between the potential wells in one subsystem is done via the parameter 
V. The parameter g still describes the particle-particle scattering in one well, but no 
scattering between the overlap of the wave functions from different wells is taken into 
account. 

The time-independent equation describing the stationary states of the complete 
system reads 

Mif) = fx%jj ( 10 ) 


with real eigenvalues /i G M because the matrix M is hermitian. Since we are also 
interested in the PT-symmetric properties of the subsystems we extend (10) to 


Mip 


Ma Me 
Mb 


i’A 



( 11 ) 


with independent eigenvalues /ij G C for both subsystems and 


Mc = 


—iy 0 
0 iy 


Mi 




( 12 ) 


for i = A, B. We can interpret (0 as two separate equations for both subsystems 
where the gain and loss is provided by the other subsystem via the matrix Mq. For 
/xa.b £ C this also allows for PT-broken states where the norm of the subsystems is 
no longer maintained, but is increased or decreased. Such solutions are therefore non¬ 
stationary states, but because the particle number of the whole system is conserved, 
there is no unlimited exponential growth or decay possible. Therefore these solutions 
describe only the onset of their growing or decaying temporal evolution. Only states 
with Ha = j^B are stationary "PT-symmetric solutions. For ha = Rb ( P^ leads to 
solutions where the gain and loss of subsystem A (represented by Im ha) is compensated 
by the loss and gain of subsystem B (Im/i^). Therefore the total particle number is 
indeed conserved. 

We can parametrize the ansatz of the wave function for this model and reduce 
the parameter count by removing a global phase. Solutions exist for different ratios of 
the probability amplitude between the two subsystems, but they may not exist over the 
whole range of the parameters. To simplify the equations we choose to restrict the norm 
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of each subsystem to one. This leads to the ansatz 

^ V’A,i \ / cos 6 *A\ 

^ V’A ,2 ^ sin 6 'Ae"^‘^A 

^ V’B,l COs6'Be+‘‘^B+i¥’rel 

t^B,2 / \ sin6'Be“'‘^®+^‘^"'=' 

with the parameters 6^ and 6b determining the distribution of the probability amplitude 
of the wave function on the two potential wells in one subsystem and the parameters (pA 
and pb describing the phase difference. The parameter p,-e\ dehnes the phase difference 
between the two subsystems. By applying additional symmetry restrictions and thereby 
reducing the parameter count even further, analytical solutions of 0 can be obtained 
and are presented in sectionAll other solutions can be gained numerically by applying 
a multidimensional root search. 


2.2. Model with a spatial resolution of the wave function 

We want to know if the basic description provided by the matrix model is sufficient 
to capture the "PT-symmetric properties and the bifurcation structure of the system 
or if a more detailed description is necessary. We do this in two steps. First we allow 
for the more detailed information of a spatially resolved wave function but retain the 
concept of isolated coupling points. The double-5 system keeps the mathematical and 
numeric intricacy at bay but still provides a spatial resolution for the wave function. 
Therefore we combine two systems with 5-potentials ([^ which describe each subsystem 
in one spatial dimension. The subsystems are then coupled at the positions of the 5- 
wells X = ±b. The depth of the potentials is controlled by . Both the depth Vq and 
the distance 2h between the wells correspond to the coupling parameter v in the matrix 
model. The coupling strength between the two subsystems is controlled by 7 and is 
only present at the two points x = ±5, i.e. the potential has no spatial extension. The 
dimensionless coupled Gross-Pitaevskii equations read 


dx'^ 




g\ipA\ + W (5(x - 5) + S{x + b)) 




+ij [5(x - h)'tl)B{h) - 5{x + h)i)B{-h)] = /iaV’a, 


dx"^ 


S'IV’bI + W - 5) + 5(x -f b)) 


fji 


-iy [5(x - b)i!A{b) - 6{x -F b)i>A{-b)] = pBi’s, (14) 

with the same physical interpretation of pa and /xb as in ( [IT| for the matrix model. 
Stationary states of the system are calculated numerically exact by integrating the wave 
functions outward from x = 0 and by imposing the appropriate boundary conditions 
on the wave functions. We require that the wave functions have to approach zero when 
X —)■ ±cx). For numerical purposes it is sufficient for the wave functions to have small 
values at X = ±Xma,r, 


x5a,b(±^ii 


0 . 


( 15 ) 
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An additional condition can be required for the norm of the wave function. In agreement 
with the normalized ansatz (13) in the matrix model we search for solutions that fulhll 

||V’a,b|P = 1- (16) 

Both wave functions are real at a: = 0. With this we enforce a global phase and the 
phase difference between the two modes at x = 0 to be (prei = 0 . 

The 10 (real) free parameters are Rep^.s, Ini/iA,s, Re'0A,B(O), Re- 0 ^ 3 ( 0 ) and 
Im0)^g(O). They are chosen such that the 10 (real) conditions, i.e. the norm (16) 
and the boundary conditions at x = ix^ax (15) are fulhlled. Note that there are no 
constraints on the /ia.b- We will see that for stationary RT-symmetric solutions the 
result is fiA = Rb ^ This is not a constraint on the root search. 


2.3. Model with a spatial resolution of both the potential well and the coupling 


We consider an additional system and remove a further restriction, viz. the point-like 
coupling approach, by duplicating the system from ([^, where the wells are formed by a 
harmonic trap and divided by a Gaussian potential barrier. This does not only provide 
us with a system with much more realistic potential wells but also allows us to extend 
the coupling of the two subsystems over the whole space. The time-independent GPEs 
of the system read 

, 02 1 . 

( ~ ^ ~ j0A + iyxe'^’^'VB = A^a^a, 

, rp, 1 . 

( ~ ^ ~ = Rsi’E- (17) 

contact ^ ^ ^ coupling 

trap 


The parameter V)p controls the height of the potential barrier between the two wells in 
one subsystem and together with the width a of the barrier it relates to the coupling 
strength v in the matrix model. Again the coupling between the two subsystems is 
controlled by a parameter labelled 7 . 

To solve this equation we use an ansatz of coupled Gaussian functions (compare 


52 53 
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5 


0 = ^ i>i,j = ^ exp {oijx^ + bijx - 1 - Cij) . (18) 

i=A,B i=A,B 

. 7 - 1,2 j - 1,2 

We use four wave functions, two for each subsystem {i = A, B) and place one in each 
well (j = 1 , 2 ). Again we place restrictions on our ansatz. We require that the norm of 
each subsystem is one, which reduces our parameter set by two. In addition we require 


ImcA,i = Pa, Imcs,i = Pe + Pveh 

IrnCAg = -pA, Imcs,2 = -PB + Pre\ (19) 

with a constant (p^ei determining the phase difference between the two modes, and again 
reducing the parameter set by two. Therefore from the 24 parameters aij,bij,Cij G C 
20 free parameters remain and must be determined such that adequate solutions are 
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found. With these constraints the ansatz is consistent with the ansatz for the matrix 
model and the system with the double-(5 potential. 

To obtain solutions of ( [Tt] ) we apply the time-dependent variational principle 
to the time-dependent GPEs 

/ o2 1 . a 

( - ^ jV’A + iya^e-^^Vn = i^V'A, 

/ o2 1 . a 

( - ^ - 9\i^+ j V’B - (20) 

contact ^ ^ coupling 

trap 

We search a parameter set for our ansatz, which minimizes the difference between the 
left-hand and right-hand side of the equation, viz. we determine the minimum of the 
functional 


54 




( 21 ) 


In this procedure '0(f) is kept constant for a given point in time and -0 = 0 is varied to 
minimize I. Since the wave function 'ip{z{t)) is not varied we require that the parameters 
z = {ajj, fcjj, do not change. A variation with respect to z leads to the equations 
of motion for the variational parameters, which follow from 


dip 

dz 


Ip — iHip ) = 0 


( 22 ) 


A more elaborate explanation of the method can be found in [^. With a numerical 
root search we can now determine those states which satisfy the 20 conditions 


0 = aij,0 = 6ij, (23) 

fii = ic* 1 = ic *2 ^ 0 = - A, 2 with i = A, B. (24) 


For PT-symmetric solutions the chemical potentials of the subsystems will fulhl ha = 
fJiB £ IK.. 


3. PT-symmetric properties and bifurcation structure of the systems 

First we will examine analytical solutions of the matrix model. The bifurcation structure 
of these solutions and their PT-symmetric properties will be discussed. Furthermore 
the differences and similarities between this four-dimensional hermitian matrix model 
and the two-dimensional matrix model with imaginary potential will be examined. 

In a next step the results obtained with the matrix model will be compared with 
the spatially extended models. Also the influence of the phase difference between the 
two modes will be investigated. 

3.1. Bifurcations structure and VT-symmetric properties of the matrix model 

To obtain analytical solutions we have to impose some constraints on the ansatz of 
the wave function of the matrix model 0- "PT-symmetric solutions must fulhl the 
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(25) 


condition Q which for this matrix model resnlts in 

with j = A, B 

and 

i’A,i = V'n,* with z = 1 , 2 . (26) 

This ensnres that the particle loss in one system is compensated by the other. These 
restrictions lead to the ansatz 

(27) 

with which we obtain an analytical expression for the chemical potentials of two VT- 
symmetric solntions 


= -|± y +7^. 


(28) 


A more detailed calculation is given in [Appendix B 


The solutions are plotted in figure and labelled si and S 2 . For different values of 
g the solutions are shifted up or down. For increasing values of 7 the difference of the 
values of the chemical potential of the two states is increased. 

■PT-broken states do not need to obey condition (25) but (26) still must be fulhlled 
since the influx and outflux between subsystem A and B must be equal. Therefore the 
ansatz for these states reads 


■0 = (cos 6 'e^‘^,sin 6 *e '“^jCos^e sin^e^^) (29) 

with fiA = R*b ■ The calculation in Appendix B yields the analytical expressions for the 
chemical potentials 

yy+w 


J‘a=7<b = -h I 2t\IP+\p'‘-P 


with P = -: 

2 


29 


(30) 


Note that the =F and ± are independent and we therefore obtain four expressions (30) for 
PT-broken states. However two of these solutions only exist in an analytically continued 
system (see figure . 

For l^fl < 2v the state S 2 passes through a pitchfork bifurcation at 


7c = 



(31) 


in which Oi and 02 are created. For 7 > 7c these two states have the same Re /ia but a 
complex conjugate Im /ta. This means that one of the states gains particles in subsystem 
A while in subsystem B it is depleted, and vice versa. The pitchfork bifurcation occurs 
at smaller values of 7 c for an increasing nonlinearity g until for \g\ = 2v the value of 7 c 
reaches zero. For values of \g\ > 2v the bifurcation between 01^2 and S 2 no longer occurs 
and the states 01^2 exist independent of S 2 for all 7 . For g < 0 the bifurcation occurs not 
with the state S 2 but with si. Thus we have shown that PT-symmetric states exist for 
the closed four-dimensional hermitian matrix model and PT-symmetry breaking can be 
observed. 
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a) 


5 = 1.5 


b) 


5 = 3.5 



■Sl ^2_Ql Q2 


Figure 3. Analytical solutions for the chemical potential of (10) are shown. The 
T^T-symmetric states are denoted by si and S 2 - T^T-broken states are labelled with oi 
and 02 . Solutions of the effective system (33) are labelled with . The states r 2 and 
r 3 only exist for jgl > 2 and therefore appear only in figure b. The coupling strength 
is set to u = 1. In a) the nonlinearity is set to 5 = 1.5 while in b) it is set to g = 3.5. 
The pitchfork bifurcation between the states 01^2 and S 2 in a) is labelled with Bp and 
occurs at 7 « 0.882. The tangent bifurcation between states r 2 and in b) is marked 
by Bt- The analytically continued solutions are plotted using lighter colours. 


Besides these states there is another class of states in the four-dimensional matrix 
model. Wave functions which fulhl the condition 

'0A,i = -iV'B.i with z = 1, 2 and V’a^, e M (32) 

lead to decoupled equations for xjjA and ipB and result in the effective two-dimensional 
model 

" 

These states effectively describe a double-well system with a real potential, where one 
potential well is lowered and the other is raised by the value of 7. As expected we hnd 
that the amplitude of the wave function in the higher well is lower than in the other. 
For 7 = 0 the system crosses over to a symmetric double-well model with no coupling 
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and therefore we can see in figure that the state ri merges with the state si and 
the state merges with the state S2- For values g > 2 v the bifurcation between the 
PT-symmetric and PT-broken states no longer exists and two new states r2 and 
emerge. Now at 7 = 0 the states ri and si as well as r2 and S2 become equal. Also r^, 
and S4 merge. For increasing 7 the states r2 and vanish in a tangent bifurcation. 


The method used to solve ( 33 ) is described in Appendix B 


We can compare the results of the four-dimensional matrix model in hgure with 
those of the two-dimensional matrix model shown in hgure [T} It is immediately clear 
that our system shows a new and richer bifurcation scenario which differs from the two- 
dimensional matrix model. While the PT-symmetric eigenvalues of the states in the 
two-dimensional system approach each other for increasing coupling strengths 7 until 
they merge in a tangent bifurcation, in our system the eigenvalues increase in distance 
for larger values of 7 and no bifurcation between the two states si^2 occurs. However, 
some generic features remain the same. In both cases the PT-symmetric state S2 with 
a real /r passes through a pitchfork bifurcation, out of which the PT-broken states with 
complex fi emerge. For both models this bifurcation moves to smaller values of 7 until, 
for a critical value of the nonlinearity g, the bifurcation vanishes and the "PT-symmetric 
and PT-broken states never coincide. 

One advantage of using a matrix model compared to systems with a more realistic 
spatially extended description is that the matrix model gives an overview over all possible 
effects in a system while remaining straightforward to calculate. Also the knowledge 
about symmetry properties and existence of states gained from the matrix model can 
help hnding states in the more complicated models, e.g. by choosing appropriate initial 
values for a root search. Since we want to concentrate our investigation on the PT- 
symmetric properties of the subsystems, we will not further investigate the states r^. 


3 . 2 . Comparison of the matrix model and the model with a spatial resolution of the 
wave function 

The results of the system with the double -5 potentials are given in hgure |^in comparison 
with those of the matrix model. To be able to compare the two models the parameters 
in the matrix model are replaced hj g ^ g/go and 7 —)■ 7/70. Also a shift A/r in the 
chemical potential is introduced. Then the parameters 70, go, v and A/r are htted to 
the results of the double -5 model. How these parameters are connected to the extended 
model can be seen in [Appendix C[ 

In contrast to the matrix model the double -5 system includes spatial properties of 
the wave functions. In hgure the wave functions for the parameters marked in hgure 
are shown. One can clearly observe the non-diherentiability of the wave functions at 
the locations x = ±5 of the 5 -potentials. It is also clearly visible that the states with 
complex chemical potential are PT-broken (see hgure [^). The two wave functions for 
the subsystems A and B fulhl the condition which ensures that the loss 

and gain in each subsystem is balanced by the gain and loss in the other subsystem 
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g = 1.5 b) g^2.0 c) ^ = 1.5, = 0.03 



0.1 0.2 0.3 0.4 0.1 0.2 0.3 0.4 0.1 0.2 0.3 0.4 


111 


Figure 4. Chemical potential /i = /iA = Mb matrix model (11) (blue dashed 

lines). The parameters of the matrix used for all three plots are = 2.75, v = 0.28 
and 7 o = 1.27. The shift of the chemical potential of the matrix model is A/r = —0.17. 
For both hgures a) and b) the phase difference (^rei was set to zero. Figure a) was 
calculated for a nonlinearity oi g = 1.5. The different states are denoted by si. 2 , 01 , 2 - 
In plot b) a nonlinearity of g = 2.0 was used. For plot c) the same nonlinearity as 
in plot a) was used but the phase difference was set to (^rei = 0.03. The figure also 
contains the results for the double-d-system (red solid lines). For the coupling of the 
two subsystems I 4 P was set to 1.0 and the d-potentials were located at 6 = ±1.1. The 
same nonlinearities as for the matrix model were used. In figure a) the parameters 
for which the wave functions are shown in figure are marked by green circles. A 
pitchfork bifurcation between the states S 2 and ai _2 is denoted by Bp. An additional 
cusp bifurcation appearing in the case i^rei is marked by Bq. 


and the PT-symmetry of the potential is maintained. Fnrthermore the wave fnnction 
of the ground state (figure [^) is much more localized in the potential wells than the 
wave function of the excited state (figure [^). 

When we compare the solutions of the matrix model with those of the model with 
the double-d potential we observe that the qualitative bifurcation structure of the states 
is the same for both models but some quantitative deviations can be seen. Before we 
continue our investigation of the cause of these differences in section 3.3 we will take a 
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ground state 


b) 


excited state 


c) 


PT-broken state 



XXX 


Figure 5. Wave functions of the double-i5 potential system for the parameter sets 
marked in figure]^, a) Wave function of the T^T-symmetric ground state, b) Wave 
function of the T^T-symmetric excited state. In c) the broken symmetry of the VE- 
broken state can be recognized. 


look at the influence of the phase difference (pre\ between the two subsystems. 

To examine the influence of the phase difference on the bifurcation scenario we 
show in figure the case in which the phase difference between the subsystems is set 
to (prei = 0.03. The pitchfork bifurcation Bp in figure turns into a cusp bifurcation 
Be- While the central ("PT-symmetric) state Si exists on both sides of the bifurcation 
point, the two outer (PT-broken) states ai,2 are created in the bifurcation of figure |^. 
In the cusp bifurcation of figure one of the outer states (depending on the sign of ^^rei) 
merges with the central state and the other outer state performs a continuous transition 
to the central state for smaller values of 7. Also the PP-symmetry of all states is broken. 
The asymmetry increases for the central state for increasing values of ^^rei- 

If we introduce the phase difference exp(i93i.ei) between the to subsystems explicitly 
into the stationary GPE (13) for the matrix model, we obtain for the subsystem A 


RAtpA,l = + ^^A,2 + sin(93i.el)7V’B,l - i COs(v3rel)7V’B,l, 

RAtpA,! = -g\'<pA,l\^tpA,2 + V^Jaa - sin((prel)7^B,2 + 1 COs((prel)7^B,2 • (34) 

'-V-' '-V-^ 

asym. pot. gain or loss 


We see that a phase difference between the two subsystems leads to different 
contributions to the real and imaginary part of the effective potential of each subsystem. 
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The real part of the effective potential can therefore become asymmetric (this not only 
depends on the phase difference ip^-ei but also on the phase value of the wave function in 
the other subsystem). 

The influence of an asymmetric double-well potential on the bifurcation structure 
has been discussed previously 55 . For an asymmetric potential there is no longer a 


pitchfork bifurcation but a tangent bifurcation. We can compare this to the well known 
normal forms of the two parameter bifurcation theory 56 . The normal form of the cusp 
bifurcation is 


0 = X = fdx) = (3 -\- ax — , (35) 

with the bifurcation parameters a and (3. In our model the role of the second parameter 
(3 is taken by the phase difference between the two subsystems. A constant ip-cei = 0 
(which is equivalent with /3 = 0) defines a line in the (prei-l parameter space. On this 
line the pitchfork bifurcation scenario emerges. 

We have seen that the phase difference between the two modes is critical to obtain 
a 'PT-symmetric system, and the breaking of this symmetry changes the bifurcation 
structure. Only for = 0 PT-symmetric states are observed. 


3.3. Comparison of the models and usefulness of the matrix model 

In the system the two modes are coupled over a spatially extended range and 
therefore the continuous change of the phase in the wave functions may play a role. 
In figure]^ we show the stationary states of the matrix model (10) in comparison with 


those of the smooth potential system (17). The parameters of the matrix model (f/oPo 


and v) and a shift of the chemical potential A/r were adjusted to the solution of the 
model 10 but remained the same for all calculations in figure with different values 
for g and (prei- In Appendix C| it is shown how the discrete matrix model can be derived 
from a continuous model. 

Again we see a pitchfork bifurcation (figure |^) in the lower state which, for 
increasing values of the nonlinearity g, moves to smaller values of 7 . The two new states 
created in this bifurcation are non-stationary (pa.b ^ I^) PT-broken states. By further 
increasing g the value of 7 at which the bifurcation occurs moves to even smaller values 
of 7 until it reaches 7 = 0 . Thus the qualitative behaviour is exactly the same as in the 
two previously investigated models. It is generic for the coupled double-well structure. 
If the phase between the two subsystems is changed to a non-zero value, the pitchfork 
bifurcation from figure]^ changes into a cusp bifurcation (compare figure]^). This is 
the same behaviour as observed in figure]^ for the double-^-potential. No change of the 
bifurcation structure or the PT-symmetric properties due to the extended coupling is 
observed. However, as can be seen in figure [^-c the agreement with the matrix model 
is nearly perfect and much better than the agreement between the matrix model and 
the model with the 5-potential wells. 

Taking a closer look at the states of the matrix model one discovers that the upper 
and lower states are symmetric with respect to —g/2 as can be seen in (28). This is no 
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9 = 0.2 *3) 0.3 5 = 0.2, (prei = 0.03 



111 


Figure 6. Comparison of the eigenvalues of the matrix model from (10) (blue dashed 
lines) with the eigenvalues of the system 0 (red solid lines), in which the BEC 
is trapped in a smooth harmonic potential separated into two wells by a Gaussian 
potential barrier. The fit parameters for the matrix model are go = 2.78, v = 0.043 
and 7 o = 0.92 and are used for all cases a)-c). The chemical potential of the matrix 
model is shifted by A/i = 2.463. The height of the Gaussian potential barrier in system 
( [TtI ) is = 0.25 with the width a = 0.5. Figures a) and c) contain the results for 
g = 0.2, while figure b) is plotted for g = 0.3. In figure c) the phase difference is 
non-zero (i^rei = 0.03). 


longer true for the models with a spatial description. To make this asymmetry visible 
we examine figure in which one state is mirrored onto the other, e.g. for one state 


/^mirror Ro R (36) 

is plotted and po is the average value of the chemical potential of both states at 7 = 0 . 
One observes that the deviation is much more pronounced in the model with h-wells 
than in the smooth potential from (17). 

In the comparison of the fit parameters Qq., v and 70 (see table [^, one parameter 
with vastly different values is evident. The coupling strength v of the two potential 
wells in the d-potential model case is approximately 6.5 times larger than in the case of 
the harmonic trap with the potential well. This means that the separation of the two 
wells is much less pronounced due to shallower wells in the case of the d-potential. This 
leads to wave functions which are not as localized as in the case of the smooth potential. 
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Figure 7. Ground state and mirrored excited state (/Tmirror = Uo — u)- The states 
are not symmetric. Figures a) and b) show the results for the Gaussian model 0 
with g = 0.2 and go = 4.854 and go = 4.2733, respectively. Figures c) and d) show 
the results of the double-d model (14) with g = 2.0 and go = —4.5 and go = —1.1, 
respectively. In the Gaussian model the hight of the potential barrier between the two 
wells in each subsystem is changed. For a) the barrier hight is Fq® = 4.0, for b) it is 
=2.5. In the case of the d-model the (real) depth of the potentials is lowered from 
Fg° = 1.0 in a) to = 2.5 in b). 


Table 1. Fit parameters of the matrix model used for the comparison with the 
spatially extended models in figures and 


Comparison with 

90 

V 

7o 

Ag 

^0 

a 

FgD b 

double-d model 

2.75 

0.28 

1.27 

-0.17 

— 

— 

1.0 1.1 

smooth potential 

2.78 

0.043 

0.92 

2.463 

2.5 

0.5 

— — 


Therefore the contribution of the overlap of the wave functions, which was negligible 
for the smooth potential, increases. The matrix model is not capable of describing the 
nonlinear interaction between wave functions of different modes. Only the nonlinear 
scattering process in the same well is taken into account. 

For further investigation one can increase the distance between the wells or deepen 
them. One might expect that the stationary states then would be in a better agreement 
with the matrix model. We compare the model with smooth potentials for different 
barrier heights (hgure Izf and hgure 0). For a lower potential barrier the asymmetry of 
the two states becomes more pronounced. The same is true for the h-model (hgure 
and hgure [^) . 
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Figure 8. Wave functions for the ground and excited states in the Gaussian model 
for different potential barriers (in a) V(p = 2.5, in b) V(p = 4.0) for a nonlinearity of 
g = 0.2. The overlap of the Gaussians at cc = 0 is much higher for the lower potential 
barrier in a) and for the excited states. 


The wave functions for the different parameter sets are shown in hgure Here the 
probability density of the ground and excited state for the smooth potential model with 
different heights for the potential barrier can be seen. One observes a higher probability 
density in the overlap region around x = 0 for the excited states. This overlap increases 
for a lower potential barrier. Thus, we can conclude that the matrix model captures 
all relevant information of the bifurcation scenario and the "PT-symmetric properties as 
long as the different potential wells are sufficiently separated. A larger overlap leads to 
quantitative changes and the loss of a mirror symmetry of pairs of values for the chemical 
potential in the (/r, 7 )-diagram, however, it does not affect the generic structure of the 
states. 


4. Summary and Outlook 

For an experimental realization of a PT-symmetric double-well potential the description 
of a physical environment which implements the gain and loss of a complex potential 
is an important prerequisite. By combining two double-well subsystems into one closed 
hermitian system we have found such a realization. 
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For the four-dimensional matrix model without a phase difference between the two 
subsystems analytical solutions for all "PT-symmetric and PT-broken states were found. 
Although the four-dimensional matrix model showed a new and different bifurcation 
scenario in comparison with the two-dimensional matrix model from 27 some generic 
features remained the same. 

The matrix model showed the same qualitative bifurcation scenario as the two 
spatially extended models. Deviations could be observed when the two wells of the 
systems were not isolated enough such that the wave functions in each well had a 
signihcant overlap between the wells. In this case the solutions from the systems with 
a spatially resolved wave function differed from those of the matrix model. A larger 
overlap leads to quantitative changes and the loss of a mirror symmetry of pairs of energy 
eigenvalues in the (/i, 7 )-diagram, however, it does not affect the generic structure of 
the states. 

The influence of the phase difference between the two subsystems was also 
examined. While the coupling strength 7 between the two subsystems took the role 
of one bifurcation parameter, the phase difference <prei took the role of another, leading 
to a two-parametric cusp bifurcation. This bifurcation degenerated for (prei = 0 into 
a pitchfork bifurcation. Only in this case PP-symmetric states could be observed 
which makes the phase difference between the subsystems critical for the PP-symmetric 
properties of the system. 

The matrix model can be investigated further. Under the assumption that the two 
wells of the system are sufficiently isolated the matrix model reduces the description of 
the system to a low number of key parameters. Therefore the analytically accessible 
matrix model of this paper could be helpful to gain more insight into the behaviour of 
coupled BECs. In particular a similar approach to realize a PP-symmetric quantum 
system via the coupling of two condensate wave functions was studied in 47 and 
revealed complicated stability properties. This system should also be representable 
in our four-mode description such that analytic expressions should be obtainable. 
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Appendix A. Time-independent solutions of the nonlinear GPE 

We will hrst consider a linear GPE in appropriate units 
d 

-\-V{x)'ip. (A.l) 

To hnd time-independent solutions one uses 

V’(t) = ^0 exp (-i/it) (A. 2 ) 

which leads to the time-indepedent equation 

/ii/>o = - Ai/>o + V {x)ipo- (A.3) 
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Solutions with a real /i are true stationary states, i.e. only the global phase is changed 
with exp(—iRe/it). By contrast states with a complex /r, in addition to the phase 
change, increase or decrease in the probability amplitude exponentially with exp(Im/it). 
For a GPE with a nonlinearity this is no longer true. If we consider 


d 

i—'ijj = -A'ljj + V{x)^p + g\^jJ\‘^^p 


(A.4) 


the previous ansatz (A. 2 ) will lead to 


IMjjo = -A'ljjQ + I/(x)V’o + 5'l^o|V^oexp(-2iIm/it). (A.5) 

Also in this case, if /i is a purely real number, the states V’ are stationary. But for states 
with a chemical potential which has an imaginary part Im/x 7 ^ 0 the interpretation 


changes. In the nonlinear case (A.2) is only a solution in the limit t —>■ 0. Therefore for 


small times the probability amplitude of these states still approximatially increases or 
decreases exponentially, but the true time evolution deviates from this linear solution 
as time increases. 

Appendix B. Analytical solutions of the matrix model 


We want to show how to calculate the analytical solution for the four dimensional matrix 
model in ( 0 . As an ansatz for PT-symmetric solutions ( [ 2 ^ is used. We obtain the 
equations 


d I — 2i(^ 

- + ve 


iye = /i. 


-1 + ne^“^ + iye^'*^ = /i. 

With the substitution x = exp(2i(p) we can transform these equations into 

(n + i 7 )x = /i+|. 

This in turn leads to 


(B.l) 


(v — iy)- = (n + iy)x 

X 


and therefore we obtain 

X = ± 


n — ly 


n + ly 


(B.2) 

(B.3) 

(B.4) 


which can be inserted into one of the equations in (B.l) and yields to the two VT- 
symmetric solutions 

/U = -| ± 


For the PT-broken solutions the ansatz (29) is used and results in 


- g cos^ 9 + V tan 9e — iye = /r, 
-g sin^ 0 + n cot 6 *e^“^ + iye^“^ = /i. 


(B. 6 ) 
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By eliminating /i and separating the equation into its real and imaginary part the 
equation system 

— 5 '(cos^ 9 — sin^ 9) + u(tan 9 — cot 9) cos 2(p = 0 

—u(tan 0 + cot 0) sin 293 — 27 cos 2(^9 =0 (B.7) 


X 


remains, which can be transformed into 

V 

sin 29 = -cos 2ip = -tan 2(p. 

9 7 

By the substitution x = exp(2i93) the quasi palindromic polynomial 

- 2Ax^ + 2x^ + 2Ax + 1 = 0 with A = -i— 

27 

is obtained. The four solutions of this polynomial are 

x=-(z±V4T^) with z = A±^A^ - 4 = with P = 

2 7 2 2g 

Note that the ± for x and P are independent and therefore lead to four solutions. By 

inserting the solutions into one of the equations in (B.6) one obtains the analytical 

expressions for the chemical potential, 


(B. 8 ) 


(B.9) 


.(B.IO) 


9 


;x = -^ 2+WP+2-P2 


T 


(B.ll) 


Note that without an analytical continuation of the equations the parameters 9 and ip 
in the ansatz of the wave functions must be real. Therefore one can see that two of the 
solutions for the chemical potential have complex 9 or ip over the whole parameter range 
and therefore are shown in lighter colours in hgure The other two solutions exist if 
the constraint 


7 > 7c = 



is fulhlled. 


For the effective matrix model in (33) with the ansatz 


-0 = (cos^, sin^) 


(B.12) 


(B.13) 


one obtains the equations 

— g cos^ 6 * — 7 + u tan 9 = jn, 

—g sin^ 0 + 7 + u cot 0 = /i, (B. 14) 

which can be transformed into the polynomial 

gy* + 4(7 + w)y^ + 4 (—7 + w)y — g = 0 (B.15) 

by eliminating fi and substituting y = e^'^. The four solutions of the polynomial can 
be obtained by any of the methods to solve polynomials of degree four. Once they are 
known ji can be calculated. 
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Appendix C. Derivation of the coefficients in the matrix model from the 
extended Gaussian model 


The matrix model ([^ can be derived as a discrete nonlinear ansatz of the extended 
model (for the derivation of a nonlinear discrete Schrodinger eqnation from the GPE 
see 46,57,|^). We rearrange the terms in (20) which resnlts in 

1 


#A = ( + -x^ + K)^e 


G ^—crx^ 


Hao 


Hi 


+ -X^ + Eo^e—^ 


Hi 


i>A+hxe >B, 

Hab 

i>B -iqxe"^""" V'A- 

Hba 


(C.l) 

(C.2) 


Also a slightly different parametrisation for the ansatz of conpled Ganssians (18) is nsed 

* = E SiJ = E = E (C,3) 

j=l,2 j=l,2 j=l,2 

with i = A, B, j = 1,2, Oij G C and Pij, qij G M. In this new ansatz the amplitnde and 
phase dij is separated from the shape fij of the wave fnnctions. It is assnmed that the 
shape is constant in time and only the amplitnde and phase changes. 

In the following paragraphs we only consider the eqnation of snbsystem A the 
calcnlation for snbsystem B can be done in the same way. We insert the new ansatz 
(G.3) into ( |G.2[ ) and obtain 

^dA,k = (Ha + Hl)dA,kfA,k + HABdB,kfB,k- (^.4) 

fc=l,2 fc=l,2 

The eqnation is multiplied with and 2 from the left and the equation is integrated 
over X. The resulting two equations can be written as a matrix equation 

(/a,i|/a,i) (/a,i|/a,2) 

(/a,2|/a,i) (/a,2|/a,2) 



=Ka =dA 

{fA,l\HA + H,\fA,i) {fA,l\HA + H,\fA,2) 
{fAMHA + H,\U,l) {fA,2\HA + H,\fA,2) 

■ V' 

=Ga 

(/a,i|77ab|/s,i) (/a,i|77ab|/b,2) 

(/Agl-f^ABl/sq) (/Agl-f^ABl/fig) 

V' 

=Gab 


dA,i 

dA,2 



(G.5) 


Gombining the equations from both subsystems results in a four-dimensional matrix 
equation 


/ Aa 0 \ \ Ga Gab W dA ^ 

Y 0 Kb j y dB j y Gba Gb j y dB j 


(G.6) 
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Table Cl. Numerical absolute values of the matrix entries for matrices Ga and Gb 


for g = 0.2 and 7 = 0.03. 

See figure]^ 


contact interaction 

external potential and kinetic energy 


Hao and Hbo 

Hi 

diagonal 

7.1214 X 10-2 

4.6221 

off-diagonal 

2.3891 X 10-^ 

8.2474 X 10-2 


We add a numerical example, which is obtained for g = 0.2 and 7 = 0.03 for the 
extended model (compare £gure|^. 

First we examine the matrices Ka and K^- 


Ka = Kb = 


1.87 

0.027+ 0.0097i 


0.027- 0.0097i 
1.87 


(C.7) 


It is obvious that the matrix has only small off-diagonal elements since the overlap 
of the wave functions of different wells is very small. Therefore the matrix K can be 
approximated by a diagonal matrix D. Then the equation is multiplied with D~^ from 
the left. 

Now we examine the matrix elements of the matrix Gab- The hrst diagonal element 
of the matrix is 


i7(/. 


A1 


xe 


—px‘ 


I/bi), 


(C.8) 


where the second term in brackets contains only structural information and can be 
integrated. It corresponds to the £t parameter 70 . We examine the numerical values of 
the matrix 

' 9.59 X 10“^+ 6.31 X IQ-^i 


Gab — 


9.40 X 10 


-5 


9.40 X 10 


-5 


9.59 X 10“^+ 6.31 X 10-£ 


- 2 \ 


.(C.9) 


It is obvious that the small overlap of the wave functions in different wells leads to very 
small off-diagonal elements, which can be neglected. 

The entries of the matrices Ga and Gb consist of terms containing the external 
potential and the kinetic energy on the one hand and terms containing the contact 
interaction on the other hand. The terms of the external potential and the kinetic 
energy in the diagonal element induce a shift of the energy (which corresponds to the 
offset A/i in the £t), therefore only the nonlinear contact interaction term remains on the 


diagonal. The contact interaction in the off-diagonal is very small (compare table Cl) 


when compared to the diagonal, and can be neglected. Therefore in the off-diagonal only 
the terms from the external potential and the kinetic energy remain. These correspond 
to the coupling parameter v in the fit. 
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